require(cjoint)
require(RColorBrewer)
Set1 <- brewer.pal(9, "Set1")
windowsFonts(Helvetica = windowsFont("Helvetica"))

attribute.label <- c("Parents", "Identity", "Activity", 
                     "Policy", "Habit", "Intention")
level.label <- list(c("Independent", "Same party", "Opponent party"), 
                    c("Does not matter", "Personal insult"), 
                    c("Never", "Rarely", "Sometimes", "Often", "Very often"), 
                    c("Does not approve", "Partially approves", "Completely approves"), 
                    c("Always abstained", "Sometimes abstained", 
                      "Sometimes deviated", "Always loyal"), 
                    c("Abstain", "Same party", "Opponent party"))

choice.formula <- formula(choice ~ a.parents + b.identity + c.activity + d.policy + e.habit + f.intention)
rating.formula <- formula(rating ~ a.parents + b.identity + c.activity + d.policy + e.habit + f.intention)

#### load dataset ####
load("entire_data.Rdata")
load("clean_data.Rdata")

# exclude Independents
entire.partisan.data <- subset(entire.data, PID != 3)
partisan.data <- subset(clean.data, PID != 3)

# number of responses
nrow(partisan.data)
round(nrow(partisan.data) / nrow(entire.partisan.data), 3)

# number of respondents
length(unique(partisan.data$respondent.id))
round(length(unique(partisan.data$respondent.id)) / 
        length(unique(entire.partisan.data$respondent.id)), 3)

#### analysis using whole partisan sample ####
AMCE.choice.total.result <- summary(cjoint::amce(choice.formula, data = partisan.data, 
                                                 respondent.id = "respondent.id"))
AMCE.rating.total.result <- summary(cjoint::amce(rating.formula, data = partisan.data, 
                                                 respondent.id = "respondent.id"))
print(AMCE.choice.total.result, digits = 1)
print(AMCE.rating.total.result, digits = 1)

png("Figure_2.png", width = 6.2, height = 3.8, units = "in", 
    pointsize = 7, res = 600, family = "Helvetica")
layout(matrix(1:2, 1, 2))
par(mar = c(2, 0, 2, 0))
plot(NULL, NULL, bty = "n", xlim = c(-1.2, 0.4), ylim = c(0, 26), 
     main = "Choice", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.4
loc.labels <- -1.07
loc.factors <- -1.2
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
  points(0, plot.location, pch = 19)
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
    points(AMCE.choice.total.result$amce$Estimate[coefficients.location], plot.location, 
           pch = 19)
    segments(AMCE.choice.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.choice.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location, 
             AMCE.choice.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.choice.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location, lwd = 1)
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.4, -0.2, 0, 0.2, 0.4))
plot(NULL, NULL, bty = "n", xlim = c(-2.4, 0.8), ylim = c(0, 26), 
     main = "Rating", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.8
loc.labels <- -2.14
loc.factors <- -2.4
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
  points(0, plot.location, pch = 19)
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
    points(AMCE.rating.total.result$amce$Estimate[coefficients.location], plot.location, 
           pch = 19)
    segments(AMCE.rating.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.rating.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location, 
             AMCE.rating.total.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.rating.total.result$amce$"Std. Err"[coefficients.location], 
             plot.location, lwd = 1)
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.8, -0.4, 0, 0.4, 0.8))
dev.off()

#### differences between PID ####
AMCE.choice.Dem.result <- summary(cjoint::amce(choice.formula, data = partisan.data, 
                                               subset = partisan.data$PID == 1, 
                                               respondent.id = "respondent.id"))
AMCE.rating.Dem.result <- summary(cjoint::amce(rating.formula, data = partisan.data, 
                                               subset = partisan.data$PID == 1, 
                                               respondent.id = "respondent.id"))
AMCE.choice.Rep.result <- summary(cjoint::amce(choice.formula, data = partisan.data, 
                                               subset = partisan.data$PID == 2, 
                                               respondent.id = "respondent.id"))
AMCE.rating.Rep.result <- summary(cjoint::amce(rating.formula, data = partisan.data, 
                                               subset = partisan.data$PID == 2, 
                                               respondent.id = "respondent.id"))

# number of responses
AMCE.choice.Dem.result$samplesize_estimates
AMCE.choice.Rep.result$samplesize_estimates
# number of respondents
AMCE.rating.Dem.result$respondents
AMCE.rating.Rep.result$respondents

png("Figure_3.png", width = 6.2, height = 3.8, units = "in", 
    pointsize = 7, res = 600, family = "Helvetica")
layout(matrix(1:2, 1, 2))
par(mar = c(2, 0, 2, 0), xpd = TRUE)
plot(NULL, NULL, bty = "n", xlim = c(-1.2, 0.4), ylim = c(0, 26), 
     main = "Choice", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.4
loc.labels <- -1.07
loc.factors <- -1.2
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
  points(0, plot.location + 0.15, pch = 19, col = 1)
  points(0, plot.location - 0.15, pch = 15, col = 1)
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
    points(AMCE.choice.Dem.result$amce$Estimate[coefficients.location], plot.location + 0.15, 
           pch = 19, col = 1)
    segments(AMCE.choice.Dem.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.choice.Dem.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, 
             AMCE.choice.Dem.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.choice.Dem.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, col = 1)
    points(AMCE.choice.Rep.result$amce$Estimate[coefficients.location], plot.location - 0.15, 
           pch = 15, col = 1)
    segments(AMCE.choice.Rep.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.choice.Rep.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15, 
             AMCE.choice.Rep.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.choice.Rep.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15, col = 1)
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.4, -0.2, 0, 0.2, 0.4))
legend(0.4, 26, legend = c("Self-identified D", "Self-identified R"), 
       pch = c(19, 15), col = c(1, 1), bty = "n", xjust = 1, yjust = 0, cex = 0.8)
plot(NULL, NULL, bty = "n", xlim = c(-2.4, 0.8), ylim = c(0, 26), 
     main = "Rating", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.8
loc.labels <- -2.14
loc.factors <- -2.4
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
  points(0, plot.location + 0.15, pch = 19, col = 1)
  points(0, plot.location - 0.15, pch = 15, col = 1)
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
    points(AMCE.rating.Dem.result$amce$Estimate[coefficients.location], plot.location + 0.15, 
           pch = 19, col = 1)
    segments(AMCE.rating.Dem.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.rating.Dem.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, 
             AMCE.rating.Dem.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.rating.Dem.result$amce$"Std. Err"[coefficients.location], 
             plot.location + 0.15, col = 1)
    points(AMCE.rating.Rep.result$amce$Estimate[coefficients.location], plot.location - 0.15, 
           pch = 15, col = 1)
    segments(AMCE.rating.Rep.result$amce$Estimate[coefficients.location] + 
               qnorm(0.025) * AMCE.rating.Rep.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15, 
             AMCE.rating.Rep.result$amce$Estimate[coefficients.location] + 
               qnorm(0.975) * AMCE.rating.Rep.result$amce$"Std. Err"[coefficients.location], 
             plot.location - 0.15, col = 1)
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.8, -0.4, 0, 0.4, 0.8))
legend(0.8, 26, legend = c("Self-identified D", "Self-identified R"), 
       pch = c(19, 15), col = c(1, 1), bty = "n", xjust = 1, yjust = 0, cex = 0.8)
dev.off()